It is well known that cell migration with its morphological dynamics is considered as a highly integrated process that is crucial for the evolution of physiological or pathological mechanisms including wound healing, immune response, embryogenesis, and cancer metastasis in biological systems.^[1–6] The investigation of cell migration from the biophysical point of view has received much attention over nearly two decades.^[6–11] It is believed that the investigations based on general physical principles will not only enrich the biophysical theory of cell interaction, but also provide a new insight and a theoretical basis for the detection^[12] and treatment^[13] of related diseases.

Cancer metastasis is one typical example of cell migration.^[14,15] Cancer cells, for example, squeeze and move through the endothelium to enter blood vessels, and eventually a secondary tumor can form in other parts of the body. Previous experiments^[16,17] have shown that the main reason for the increased mobility of migratory cancer cells is the difference in mechanical properties between cancer cells and normal cells, i.e., cancer cells show stronger intercellular adhesion and protruding activity. Recent experiments by Lee et al.^[16] showed that the Young modulus of human breast epithelial cells (MCF10 A) is approximately three times larger than that of metastatic human breast carcinoma cells (MDA-MB-231). They observed that the mobility of a cancer cell embedded in a confluent monolayer of normal cells is much greater than that in a confluent monolayer of cancer cells. From a mechanical point of view, the cancer cell is stiffer than the normal cell, hence cancer cells move faster than normal cells. Therefore, understanding the cell migration behavior with shape change in the context of their mechanical properties is an important biophysical problem. However, it is hard to experimentally isolate the role of elasticity mismatch on cell migration in experiments.

With the rapid development of computational mathematics and computer technology, numerical simulation has gradually become a scientific research method which can be parallel to experimental and theoretical research. Dynamic changes of cells in shape are the obvious difference between cell and particle migration. Cell migration dynamics with morphological changes can be considered as an out-of-equilibrium system, which can be mathematically described as a free boundary problem. Many numerical attempts^[11,18–25] have been made to cell migration, such as isotropic particles, simple deformable particles, vertex models, subcellular element models, and cellular Potts models.

Due to the ability to avoid the complex boundary tracking algorithm, the phase-field method over decades has gradually become a well-accepted technique to solve free boundary problems by introducing an order parameter “phase-field”.^[26–34] Recently, the phase-field method has been extended to cell migration and morphological changes in biological fields. Shao et al.^[4,35] developed a phase-field model incorporating the membrane bending force and surface tension to investigate single cell morphodynamics. The predicted steady-state cell shapes with a wide range of aspect ratios depend on system parameters, agreeing quantitatively with experiments. Based on a Ginzburg–Landau free energy formulation, Najem and Grant constructed a series of phase-field models for neutrophil morphodynamics,^[36] neural cell chemotropism,^[37] an immune response,^[38] and collective cell migration,^[39] showing that the phase-field method is a powerful and general approach to simulate cell migration. By using a developed phase-field model for a crawling cell, Löber et al.^[40] simulated the movement of crawling cell on soft substrates, showing complex cell movement modes. Using a phase-field to account for the bending elasticity of the membrane, Lázaro et al.^[41,42] investigated the shapes of red blood cells under flow in highly confined micro-channels; the results showed that the relevant aspects of cell elasticity contribute to the shape of blood under high confinement, and the increase in the flow velocity results in a sharper cell distribution in the channel. Very recently, Palmieri et al.^[17] proposed a multi-phase-field model for monolayer of motile cells comprising normal and cancer cells. This model focused on the role of the elastic properties of cancer cells and surrounding cells on the self-migration of cancer cells, and simulated the mobility of a cancer cell embedded in a confluent monolayer of normal cells. Their results showed that the elasticity mismatch alone is sufficient to increase the motility of the cancer cell, and the increased motility and the amplitude and frequency of the bursts agree quantitatively with experiments. However, the mechanisms of these interactions among cells are not clearly understood yet, and cell migration behaviors with different mechanical properties in a confluent monolayer have not been addressed in detail. A more precise understanding of the elastic mismatch on the migration dynamics of multiple cells in a confluent monolayer is highly desirable. Therefore, in this paper, the phase-field model by Palmieri et al.^[17] is used to exhibit the effect of elasticity mismatch on the migration dynamics of multiple cells in a confluent monolayer. Through a large number of numerical simulations, the effect of cell elasticity mismatch and the overall migration speed on the cell shapes has been investigated. The remainder of this work is organized as follows. In Section 2, we briefly describe the phase-field model for multi-cell migration and the parameter selection. We present our simulated results and discussions in Section 3. Section 4 contains conclusions from this work.

We used the phase-field model for cell migration dynamics in monolayer by Palmieri et al.^[17] to investigate the effect of elasticity mismatch on the cancer cell migration in a confluent monolayer of normal cells. Compared with other models for cell migration, no need to explicitly track the cell boundary and the ability to describe of the extremely large deformation are the advantages of the phase-field model. Moreover, the phase-field model^[17] can be considered as a more general model for cell migration, because (I) it is easy to extend the phase-field model to more complex cases (such as three-dimensional system, intracellular dynamics, etc.); (II) the dynamics of the phase-field model is self-consistently derived from non-equilibrium thermodynamic principles; (III) the phase-field model focuses on the single cell behavior, which is important to describe the interactions among the normal cells and the tagged cells. The basic idea of the phase-field model is that a continuous variable, “phase-field” φ_i(r;t), is introduced to distinguish the interior of the cell or its surrounding, where the subscript i denotes the cell number. The phase-field variable φ_i (r;t) takes on φ_i(r;t) = 1 in the interior of cell i, and φ_i(r;t) = 0 represents the region outside the cell, and it varies rapidly from 1 to 0 in a numerically resolvable diffusion interface region, within which the cell boundary can be considered as the contour of φ_i(r;t) = 0.5. Therefore, the boundary of the cell i can be implicitly tracked by φ_i(r;t). The migration of each cell in the monolayer in an n-cell system can be described by the dynamic equation of phase-field variable, yielding 1 where Γ_i is the dynamic mobility, F is the total free energy which can be divided into two parts, F = F_i + F_int, where F_i is the free energy for individual cell, and F_int is the free energy for the interaction between cells. In a two-dimensional (2D) system, F_i and F_int can be written as 2 3 where λ is related to the width of the diffuse cell boundary, γ_i is defined to describe the elasticity of the cell, R is the radius of the non-interacting cell, and κ controls the energy cost resulting from the overlap of the cell N and cell M. Equation (1) describes the dynamic migration of each cell in a monolayer, and hence the translational velocity of the cell i is very important for this work. The translational velocity of the cell i can be divided into two parts v_i = v_i,a + v_i,ia, where v_i,a and v_i,ia are active part and inactive part of the velocity, respectively. The active part is owing to the internal process requiring energy cost, while the inactive part is due to the interactive force by the neighboring cells, yielding 4 where ε determines the relative strength of the two parts to the velocity. The dimensional form of the cell dynamic equation can be written as 5 6 Here, the prime symbol (′) denotes the dimensionless variable, λ′ = λ/Δx, , , κ′ = λ′²Δx²κ/(30γ₀), μ′ = πR′²Δx⁴μ/γ₀, t′ = 2Γ_it/Δx², and ε′ = Γ_iε, where Δx is the grid spacing, and γ₀ is the cell stiffness of normal cell. For details of this phase-field model, interested readers can refer to Ref. [17].

In the present work, we aim at investigating the effect of cell elastic mismatch on the migration of multiple cells in monolayer. To study this, two types of cells with different cell stiffness are considered, i.e., one tagged cell and n − 1 surrounded normal cells, as shown in Fig. 1. In this model, γ_normal and γ_tag are defined to describe the elasticity of the normal cell and the tagged cell, respectively. For normal cells, γ_normal is fixed at 1; while for the tagged cell, γ_tag varies from 0.1 to 1.9 for each simulation. Note that γ_tag / γ_normal < 1 means that the tagged cell is softer than the surrounded normal cells; on the contrary, γ_tag / γ_normal > 1 means that the tagged cell is harder. Different types of cells have different elasticity ratios, e.g., γ_tag / γ_normal = 0.35 is determined by the in-situ observation in experiments in a cancer–normal system.^[16] All other parameters are assumed constants^[17] and listed as follows: λ′ = 7, R′ = 12, κ′ = 60, μ′ = 40, and ε′ = 1.5 × 10³. |v_a| varies from 0.1 to 0.3, where |v_i,a| is replaced by |v_a| for simplicity because all |v_i,a|′s are the same in all simulations. The simulations are performed in a square simulation domain with 200 × 200 mesh grids, and periodic boundary conditions are imposed. The simulations are initialized by randomly setting n cells in the simulation domain, including one tagged cell and n − 1 normal cells with radius R′. The time-dependent dynamic equations (Eqs. (5) and (6)) are solved by using the explicit finite difference method on a fixed grid and time stepped with a simple Euler scheme. Moreover, it should be noted that the simulation results can be compared with experiments by comparing the velocity bursts in simulations and experiments, although a statistically meaningful comparison between predicted results and experimental observations is difficult to be done because of limited experimental observations, as discussed in Ref. [17].

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) A schematic of migration dynamics of multiple cells in a confluent monolayer, where a tagged cell (red one) is embedded into an array of normal cells (white ones) and periodic boundary conditions are imposed.

In the initial stage of simulations, we randomly set n = 65 cells that are nearly circular in the simulation domain with radius R′ = 12, which means that the cell packing fraction is fixed at 0.735. Following the method given by Ref. [17], we solve Eq. (5) without the active part of the velocity (|v_a| = 0) at the initial stage because many overlapping cells are far from equilibrium. This initial stage is the so-called “aging” period. Figure 2 shows the snapshots of motile multiple cells including the tagged cell and other normal cells for |v_a| = 0.3 (a) in the soft-in-normal case (γ_tag / γ_normal = 0.4), (b) in the nearly all-normal case (γ_tag / γ_normal = 1.1), and (c) in hard-in-normal case (γ_tag / γ_normal = 1.9), where each cell is described by the phase-field variable and the green one is the tagged cell. Here, the normal cells are taken as the same, and hence we artificially add a small amount (0.1) to the phase-field variable for the tagged cell to distinguish the tagged cell from others. In the nearly all-normal case (γ_tag / γ_normal = 1.1), since γ_tag is very close to γ_normal, the deformation and migration behavior of the tagged cell are nearly similar to those of normal cells, as shown in Fig. 2(b). It is obvious that the tagged cell squeezes and migrates with large deformation among normal cells in the soft-in-normal case because the tagged cell is softer than the surrounded normal cells for γ_tag / γ_normal = 0.4, as shown in Fig. 2(a). On the contrary, figure 2(c) shows that the tagged cell moves slowly with small deformation and normal cells migrate rapidly around it for γ_tag / γ_normal = 1.9. Therefore, the cell elasticity mismatch indeed plays an important role in multi-cell dynamics.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Snapshots of phase-fields for motile multiple cells including one tagged cell and other normal cells for |va| = 0.3 at different elasticity ratios γtag / γnormal (a) 0.4, (b) 1.1, and (c) 1.9, where the time interval between each snapshot is 50 and a small value (0.1) is added into the phase-field variable for tagged cell to distinguish the tagged cell from the normal cell.

In these simulations, the volume (area in 2D) of the cell is assumed to be constant. In order to characterize the deformation and distortion of the cells, the dimensionless length of the cell perimeter is used here, and it can be calculated by 7 It can be seen from the definition of Eq. (7) that the shapes of the cells become perfectly circular when while the deformation becomes large with the increase of . To compare the deformation of tagged cell with those of normal cells, we plot the variations of perimeters of the tagged cell and arbitrary four chosen normal cells with time for the three cases (γ_tag / γ_normal = 0.4, γ_tag / γ_normal = 1.1, and γ_tag / γ_normal = 1.9), as shown in Fig. 3. In the soft-in-normal case (γ_tag / γ_normal = 0.4), the difference between tagged cell and normal cells is very intuitive. The tagged cell undergoes very large deformation compared with the chosen four normal cells: the maximum value of the perimeter of the tagged cell exceeds 2.0, i.e., 200% of the initial circular shape, while those of the normal cells do not exceed 1.4, as shown in Fig. 3(a). In the nearly all-normal case (γ_tag / γ_normal = 1.1), it is difficult to distinguish the tagged cell from normal cells. As shown in Fig. 3(b), the difference of perimeters among tagged cell and chosen normal cells is very small as the cell elasticity mismatch is relatively small. Although the tagged cell is a little harder than the normal cells, the normal cells have relatively large deformation and the perimeter can exceed 1.4 in this case while the perimeter of the tagged cell just exceeds 1.3. As for the hard-in-normal case (γ_tag / γ_normal = 1.9), except for the initial stage, the perimeter of the tagged cell keeps at a very low value, while the maximum value of the perimeters of the chosen normal cells can further increase and exceed 1.6, as shown in Fig. 3(c).

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) The dimensionless perimeter L′ of the tagged cell and four normal cells as a function of dimensionless time for (a) γtag / γnormal = 0.4, (b) γtag / γnormal = 1.1, and (c) γtag / γnormal = 1.9, where the solid line corresponds to the tagged cell and dashed lines correspond to normal cells.

Figure 4(a) shows that the evolution of the dimensionless instantaneous velocity |v| and the dimensionless perimeter L′ of the tagged cell for γ_tag / γ_normal = 0.4. The curve of the instantaneous velocity versus time shows oscillations with some short speed bursts. The instantaneous velocity |v| exceeds 2, which is about ten times of |v_a|. This indicates that the deformation can provide a large velocity of the cell. To clearly investigate the relationship between the instantaneous velocity and the perimeter, a magnification of Fig. 4(a) from t = 1200 to t = 1400 is plotted in Fig. 4(b). It can be observed from this magnification that the instantaneous velocity firstly decreases with the increase of the perimeter, and then the short speed burst occurs simultaneously with a significant decrease of the perimeter. As the elasticity energy releases and is converted to kinetic energy, the cell quickly becomes relaxed and its shape becomes more circular. Therefore, the bursts of the velocity always fall behind the peak of the perimeter. Figure 5 illustrates the correlation between the dimensionless instantaneous acceleration d |v|/dt and cell perimeter change ΔL′ for γ_tag / γ_normal = 0.4 and for γ_tag / γ_normal = 1.9. The instantaneous acceleration d |v|/dt versus cell perimeter change in long simulations is tracked. Generally, the larger the cell perimeter change, the faster the instantaneous acceleration. It is found that the distributions of the values of the instantaneous acceleration according to the cell perimeter change for the two γ_tag / γ_normal values are very close, showing that the shape change of the tagged cell results in the accelerated motion of the tagged cell. This is irrelevant to the elasticity mismatch.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) (a) The dimensionless perimeter L′ for the tagged cell and dimensionless instantaneous velocity |v| of the tagged cell as a function of dimensionless time for γtag / γnormal = 0.4, where the red line corresponds to the dimensionless perimeter and black line corresponds to the dimensionless instantaneous velocity; (b) magnification of (a) from t = 1200 to t = 1400.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) The dimensionless instantaneous acceleration d |v|/dt versus cell perimeter change ΔL′ within a long time for γtag / γnormal = 0.4 and for γtag / γnormal = 1.9.

Then, we analyze the variation of the shape change of cells with the variation of γ_tag / γ_normal for different |v_a| values. To explore the more universally applicable law of the cell perimeter with the elastic mismatch, γ_tag / γ_normal is chosen in a wide range from 0.1 to 1.9, i.e., from soft tagged cell to hard tagged cell. We average the shape variables (the cell perimeter) for tagged cell and normal cells under various conditions and obtain the average cell perimeter for tagged cell and normal cells. In order to eliminate the random errors caused by different cell shape variables in each simulation, we use the cell perimeter ratio, , to characterize the effects of elasticity mismatch on cell shapes with the variation of γ_tag / γ_normal for three different |v_a| values, and the results are shown in Fig. 6. For the soft tagged cell, the shape change ratio is greater than 1; when γ_tag / γ_normal > 1 (hard tagged cell), is lower than 1. It can be seen that is a decreasing function of γ_tag / γ_normal, which can be described by a simple power law function for three values of |v_a|, implying the active velocity of the cell is not an important factor for cell shape change, at least in this range.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) The average shape change of tagged cell over the average shape change of normal cell as a function of γtag / γnormal for three different |va| values.

In summary, we use the phase-field model for cell migration to simulate the migration dynamics of multiple cells in a confluent monolayer with elastic mismatch, where one tagged cell is embedded by a number of normal cells and both types of cells have the same properties except elasticity. Our results show that the cell elasticity mismatch plays an important role in multi-cell dynamics. Larger elasticity mismatch results in larger shape change difference between the tagged cell and the normal cells. We find that the bursts of the velocity always fall behind the peak of the perimeter, and the shape change of the tagged cell results in the accelerated motion of the tagged cell in the whole process. Moreover, the variation of the shape change of cells with γ_tag / γ_normal for different |v_a| values is analyzed. We find that is a decreasing function of γ_tag / γ_normal, which can be described by a simple power law function. In future work, by considering the volume change dynamics of cell, this phase-field model can be extended to the case of the multiple tagged cells-multiple normal cells migration dynamics for investigation of tumor growth and invasion in a larger scale.^[43]